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I extend the scope of the density matrix renormalization group technique developed by White to 
the calculation of dynamical correlation functions. As an application and performance evaluation I 
calculate the spin dynamics of the ID Heisenberg chain. 



The density matrix renormalization group method 
(DMRG) as developed recently by S. White @] is a pow- 
erful algorithm for calculating ground state energies and 
static properties of low dimensional systems. This tech- 
nique leads to highly accurate results for much larger sys- 
tems than those which can be solved by straightforward 
exact diagonalization. The method has been applied suc- 
cessfully to several problems such as the Haldane gap of 
spin-1 chains the one-dimensional Kondo insulator ||] 
and the two-chain Hubbard model Q . 

An effective way of extending the basic ideas of this 
method to the calculation of dynamical quantities was 
lacking, mainly due to the fact that it is performed in 
real space and it is not possible to fix the momentum as 
a quantum number ||. It also involves a strong trunca- 
tion of the Hilbert space and therefore much information 
of the excited states is lost. 

In this paper I present a way to calculate the dynam- 
ical properties using the DMRG method. As an appli- 
cation I calculate the spin dynamics of the ID isotropic 
Heisenberg model. The dynamics of this model has been 
studied extensively [PHlij and therefore presents a good 
background for comparison. 

The DMRG allows for a systematic truncation of the 
Hilbert space by keeping the most probable states in de- 
scribing a wave function (e.g. the ground state) of a larger 
system, instead of the lowest energy states usualy kept 
in previous real space renormalization techniques. The 
method is very well described in Ref. [Q but I shall sum- 
marize it so as to unify notations. A general iteration 
of the method consists of: i)The effective Hamiltonian is 
defined for the superblock l+2+l'+2' (a block is a col- 
lection of sites), where the blocks 1 and 1' come from 
previous iterations and blocks 2 and 2' are new added 
ones. It is diagonalized to obtain the ground state |^>o) 
(other states could be also kept: they are called target 
states), ii) The density matrix pa> = J^. tpo^jtpo^'j is 
constructed, where Vo,ij = (* ® j'lV'o); the states \i) (|j)) 
belonging to the Hilbert space of blocks 1 and 2 (1' and 
2'). The eigenstates of p with the highest eigenvalues 
(equivalent to the most probable states of blocks 1+2 in 



the ground state or in the chosen target state of the su- 
perblock) are kept up to a certain cutoff, keeping a total 
of m states per block, iii) These states form a new re- 
duced basis to which all the operators have to be changed 
and the block 1+2 is renamed as block 1. iv) A new 
block 2 is added (one site in our case) and the new su- 
perblock (l+2+l'+2') is formed as the direct product of 
the states of all the blocks (the blocks 1' and 2' are iden- 
tical to blocks 1 and 2 respectively). When more than 
one target state is used, i.e more than one state is wished 
to be well described, the density matrix is defined as: 



Pi 



(i) 



where pi defines the probability of finding the system in 
the target state \4>i) (not necesseraly eigenstates of the 
Hamiltonian). 

I want to calculate the following dynamical correlation 
function at T = 0: 



C A {t-t') = (^\A\t)A{t')\^), 



(2) 



where A^ is the Hermitean conjugate of the operator A, 
A(t) is the Heisenberg representation of A, and |-0o) is 
the ground state of the system. Its Fourier transform is: 

C A (w) = ]T \(t/j n \A\i,o)\ 2 S(oj - (E n - E )), (3) 



where the summation is taken over all the eigenstates 
\tjj n ) of the Hamiltonian H with energy E n and Eq is the 
ground state energy. 

Defining the Green's function 

G A {z) = (MA\z-H)- 1 A\i> ), (4) 

the correlation function CU(w) can be obtained as 

CUM = -- lim Im Ga(w + ii] + E Q ). (5) 

The function Ga can be written in the form of a con- 
tinued fraction: 



1 



G A (z) 



(G) 



z — a,Q 



The coefficients a n and b n can be obtained using the fol- 
lowing; recursion equations 



\fn+i)=H\f n )-a n \f n )-b 2 n \f n - 1 ) 



(7) 



where 



l/o) = A|Vo) 
a n ={fnW\fn)/(fn\fn), 
b n = (/n|/n)/(/n-l|/n-l)j 



60 = o 



(8) 



An alternative way for calculating the spectra is by 
means of the Liouvillian representation of the recursion 
method presented above 15 13] . This method leads to 
quasi size-independent coefficients. In the example given 
below, it has been seen that the results are the same as 
with the Hamiltonian representation using Eqs. (7) and 
(8) f§. 

For finite systems the Green's function Ga(z) has a 
finite number of poles so only a certain number of coeffi- 
cients a n and b n have to be calculated. The DMRG tech- 
nique presents a good framework to calculate such quan- 
tities. With it, the ground state, Hamiltonian and the 
operator A required for the evaluation of Ca(w) are ob- 
tained. An important feature is that the reduced Hilbcrt 
space should also describe with great precision the rele- 
vant excited states \ip n )- This is achieved by choosing the 
appropriate target states. For most systems it is enough 
to consider as target states the gound state Vo) and the 
first few |/„) with n — 0,1... and |/o) = A\i/)o) as de- 
scribed above. In doing so, states in the reduced Hilbert 
space relevant to the excitated states connected to the 
ground state via the operator of interest A are included. 
The fact that |/o) is an excellent trial state, in particular, 
for the lowest triplet excitations of the two-dimensional 
antiferromagnet was shown in Ref. fl7|] . 

Another straightforward election of target states is to 
take the excited eigenstates and the ground state. But 
this is possible only when the quantum number of the 
above mentioned states can be fixed. Otherwise, when di- 
agonalizing the Hamiltonian, the lowest lying excitations 
of the whole space are obtained and not those of the ap- 
propriate symmetry sector. With the DMRG technique 
we can fix the total and parity but not, for example, 
the momentum q of the system. In the example given 
below, I attempted to obtain the first excited states for a 
given q by using |/ ) = S*\ip ) (where = £\ e iqR i S 7 -) 
as a trial state for the diagonalization procedure in step 
(i) above. Due to the fact that there is a small but non- 
zero overlap between |/o) and the algorithm does 
not remain in the symmetry sector with a given q. Be- 
cause of this, I found it convenient to use |^o) an d \fn) 
with n — 0, 1.. as target states. 



Of course, if the number m of states kept per block is 
fixed, the more target states considered, the less precisely 
each one of them are described. An optimal number of 
target states and m has to be found for each case. Due 
to this reduction, the algorithm can be applied up to 
certain lenghts, depending on the states involved. For 
longer chains, the higher energy excitations will become 
inaccurate. Proper sum rules have to be calculated to 
determine the errors in each case. 

As an application of the method I calculate 

S"(q,u) = £ K^l^lV'o)! 2 6(w - (En - E )), (9) 

n 

for the ID isotropic Heisenberg model. 

As I already mentioned, the spin dynamics of this 
model has been extensively studied. The lowest ex- 
cited states in the thermodynammic limit are the famous 
des Cloiseaux-Pearson (dCP) triplets jy], having total 
spin S T = 1. The dispersion of this spin- wave branch is: 



-7r|sin(?)l 



(10) 



Above this lower boundary there exists a two-parameter 
continuum of excited triplet states that have been calcu- 
lated using the Bethe ansatz approach fisfl with an upper 
boundary given by 



= j7r|sin(q/2) 



(11) 



It has been shown Q , however, that there are excitations 
above this upper boundary due to higher order scatter- 
ing processes, with a weight that is at least one order of 
magnitude lower than the spin-wave continuum. Based 
on selection rules, Bethe ansatz calculations and numer- 
ical diagonalization of small clusters, Miiller et al. || 
proposed the following approximate expression for the 
out-of-plane dynamical structure factor: 



.4 



=e(w-^)e«-«) (12) 



S zz (q,u) = 



where A is a constant and Q(x) a cutoff step function 
that was considered so that the sum-rules are satisfied. 
A similar expression for an exactly solvable model (the 
Haldane-Shastry model) has been obtained jlj| . The low 
energy properties of this model and those of the nearest 
neighbour Heisenberg model we are considering belong 
to the same universality class. 

In the following I will present the numerical results. 
The values for TV = 20 sites (without reduction of the 
Hilbert space) and exact calculations using the Lanczos 
technique and exploiting all the symmetries coincide ex- 
actly. For larger systems I used m = 200 states per block 
and periodic boundary conditions. 

In Fig. 1 I show the spectrum for various systems 
lenghts and q — ir and q = tt/2. The delta peaks of 



2 



Eq. (|9|) are broadened by a Lorentzian for visualizing pur- 
poses. For this case it was enough to take 3 target states, 
i. e. \ipa), |/o) = S*\tpo) and I also plot the an- 

alytical expression (|T2|). There is good agreement up to 
TV ~ 40 with the envelope of our data. For larger values 
the peaks at w/J ~ 2 acquire high weight, which grows 
with TV. The second peak seems to be somewhat shifted, 
also having a higher weight. Due to the truncation of the 
Hilbert space the spectrum of H is also reduced. I notice 
that if we consider only the first (~ 10) coefficients a n 
and b n , the spectrum at low energies remains essentially 
unchanged. Minor differences arise at uo/ J ~ 2. This is 
another indication that only the first |/„) are relevant for 
the low energy dynamical properties for finite systems. 

In the inset of Fig. 1 the spectrum for q — n/2 and 
TV = 28 is shown. For this case I considered 5 target 
states i. e. \ip ), |/ ) = S* /2 \ip ), \f n ) n = 1,3 and 
rn = 200. Here, and for all the cases considered, I have 
verified that the results are very weakly dependent on the 
weights pi of the target states, as long as the appropriate 
target states are chosen. For lenghts where this value of q 
is not defined I took the nearest value. I found that with 
these parameters the spectrum starts developing spurious 
peaks for larger systems. The coefficients also present a 
larger dispersion with TV than for the q = n case. To be 
able to go further, one should consider larger m values. 

The q = n/2 case is the most unfavourable one because 
it involves high energy excitations. Although I have in- 
cluded some of these states as target states, the reduced 
Hilbert space related to |^o) and S z , 2 |-0o) have very small 
overlap and many states are needed to describe correctly 
both target states. For q — tt, instead, the overlap of 
Hilbert spaces is very high and the target states and low 
energy excitations are better described. The difference 
between the q = n/2 and q = tt cases can be seen by 
comparing the ground state energy, it being more precise 
in the latter case by a factor of 3 in the relative error for 
m = 200 and TV = 28 (the relative error of the ground 
state for the q = tt case is 10~ 6 ). 

Even though I am including states with a given mo- 
mentum as target states, due to the particular real-space 
construction of the reduced Hilbert space, this transla- 
tional symmetry is not fulfilled and the momentum is 
not fixed. To check how the reduction on the Hilbert 
space influences the momentum q of the target state 
|/o) = S*\ipo), I calculated the expectation values 

(ih\S^SI\ih) (13) 

for all q' . If the momenta of the states were well de- 
fined, this value is proportional to 6 q - q i if q ^ 0. For 
q = 0, y\ S z = 0. In Fig. 2a) I show the expectation 
values ( |l3|) for q = n/2 (using S^, 2 \tpo) as one of the 
target states) and different lenghts. I see that, as the 
system becomes larger (higher reduction of the Hilbert 
space), the q value becomes less defined, presenting a 



wider distribution. The figure shows a marked oscilla- 
tion in the expectation value. This is due to the fact 
that the system is built from two identical blocks and 
that KV'ol'SfflV'o)! is small but non-zero (~ 10~ 3 for the 
largest system). These should disappear when using the 
finite-size method pL The momentum distribution for 
q = tt is shown in Fig. 2b) in a semilogarithmic scale. In 
this figure I have shifted the values by .003 so as to have 
well-defined logarithms. I have also neglected the points 
where Eq. ( |l3| ) is zero, mentioned above (i. e. between 
any two successive values in the figure there is a q' that 
leads to a zero expectation value) . I can see here that the 
momentum is better defined, even for much larger sys- 
tems, but, as expected, more weight on other q' values 
arise for larger TV. I also calculated Eq. ( |l3| ) for TV = 28 
and q = n/2 but using S%\tpa) as a target state. I find 
a very poorly defined momentum centered at q' = n/2. 
This is expected since the reduced Hilbert space targeted 
q = n states (in addition to the ground state). 

In Fig. 3, I show the dispersion curve for 28 sites as 
compared to the exact dCP dispersion (Eq. (|To|)). The 
difference in the values at higher values of q is due to 
finite-size effects. These results are in very good agree- 
ment with those of Ref. [^0|. The inset shows the first 
excitation energy w 1 for q = n as a function of 1/TV. In 
the thermodynamic limit this value must go to zero. The 
upwards curvature for the largest systems is due to the 
approximation of the method. 

I find excellent agreement in excitation energies and 
weights for all values of q with exact results for TV = 24 
sites 

As a check of the approximation I calculated the sum 
rule 

i r°° r 27T i 
^ jf J S"(q,w) = (M^ =0 ) 2 hV = 4 (14) 

for TV = 28, 5 target states and m = 200. I obtain a rela- 
tive error of 0.86%. I have also found that the expression 
for the static structure factor S zz (q) given in Ref. |h| 
i.e. S zz (q) = — 1/4 ln(| 1 — q/n\), fits very accurately our 
data for TV = 28 (details will be given elsewhere). 

To conclude, I have developed a method to calcu- 
late dynamical correlation functions precisely using the 
DMRG technique to evaluate the coefficients of the con- 
tinuous fraction representation of the Green's function. 
I show that even by considering a 0.1% of the total 
Hilbert space (for TV = 28 only ~ 40000 states are kept) 
a reasonable description of the low energy excitations 
is obtained. I also show that it is possible to obtain 
states with well defined momenta if the appropriate tar- 
get states are used. The numerical computation has 
been performed on a workstation. A better performance 
(more accuratetly described excitations, larger systems) 
can surely be obtained by supercomputers, where, due 
to a larger memory space, a larger reduced Hilbert space 
can be considered. 
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FIG. 1. Spectral densities for q = n and N = 28 (contin- 
uous line), N — 40 (dotted line) and N = 60 (dashed line). 
The analytical expression (12) is shown with a dashed-dotted 
line. Inset: Spectral density for q — tv/2 for N = 28 (n = .05). 



FIG. 2. Momentum weights (Eq. Ql_^) ) of a target state 
with a) q — tt/2 and N — 24 (circles), N — 28 (squares) and 
N = 36 (diamonds); b) q = n for N = 28 (circles), N = 44 
(squares), N = 60 (diamonds) and N — 72 (triangles). The 
dotted lines are a guide for the eye. 



FIG. 3. Dispersion relation of the lowest energy excita- 
tion. The full line is the spin wave curve ui l q (Eq. (jnj)). 
The dots are our data for N — 28 and m = 200. In- 
set: First excitation energy as a function of 1/JV for 
N = 14, 20, 28, 32, 36, 40, 44, 52, 60 and 72. The dotted line 
is a guide to the eye. 
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Fig. 2b 
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